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A condensate in an optical lattice, prepared in the ground state of the superfluid regime, is 
stimulated first by suddenly increasing the optical lattice amplitude and then, after a waiting time, 
by abruptly decreasing this amplitude to its initial value. Thus the system is first taken to the 
Mott regime and then back to the initial superfluid regime. We show that, as a consequence of this 
nonadiabatic process, the system falls into a configuration far from equilibrium whose superfluid 
order parameter is described in terms of a particular superposition of Glauber coherent states that 
we derive. We also show that the classical equations of motion describing the time evolution of this 
system are inequivalent to the standard discrete nonlinear Schrodinger equations. By numerically 
integrating such equations with several initial conditions, we show that the system loses coherence, 
becoming insulating. 
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1. INTRODUCTION 

Nowadays Bose-Einstein condensates (BECs) repre- 
sent one of more powerful and versatile testing grounds 
for low-energy modern physics in which experimental 
tests on quantum computation [![, many-body physics 
0, superfluidity the Josephson junction effect @, 
atom optics Q , and quantum phase transitions Q can be 
performed. Condensates can be put into interaction with 
each other or manipulated by means of optical lattices 
(OLs), which are periodic trapping potentials generated 
by standing laser waves. By raising the amplitude of the 
OL, a condensate loaded therein is fragmented into an 
array of interacting condensates. By adjusting the laser 
amplitudes, the system is taken into different regimes. 
The superfluid regime is obtained with weak optical po- 
tentials (OPs), where the kinetic energy dominates over 
the interacting one, and the atoms hop from one well to 
another. The opposite -quantum- regime, is generated 
by strong OPs that suppress the tunneling of the atoms 
between the wells. 

In this paper, we consider a one-dimensional gas. This 
is prepared by use of a transverse harmonic confining 
potential that tightly confines the atoms so that their 
motion, in the transverse direction, is limited to the zero 
point. Along the longitudinal direction, a further har- 
monic potential weakly confines the atoms, and an OP is 
switched on. For large enough laser amplitudes, the con- 
densate splits into components tightly confined at the 
minima of the effective potential. In what follows, we 
imagine abruptly adjusting the amplitude of the longi- 
tudinal laser at two instants, in order to first take the 
system from the superfluid regime to the quantum one 
and then to take it back to the superfluid regime. In this 



way the system is taken to a nonequilibrium state (NES) 
that is described in terms of a particular superposition 
of Glauber coherent states (CSs). The latter combina- 
tion is derived in the following. Although the system is 
taken back to a weak OP regime, as in the superfluid 
case, this NES follows a nonclassical dynamics. In fact, 
we show that the equations of motion for the system's 
order parameter are inequivalent to the discrete nonlin- 
ear Schrodinger equation. The present study is in some 
sense complementary to former work [9j, lldj , it continues 



the work in 
as in Refs. 




and it is motived by recent experiments 
12l . [l3| , where BECs are manipulated in 
an OL. Furthemore, the nonadiabatic procedure we de- 
scribe, can be straightforwardly realized in a real exper- 
iment, and then the dynamics of the NES, a superposi- 
tion of the canonical CSs can be directly observed, e.g., 
by displacing the condensates with respect to the har- 
monic trap center, as in the experiment of Ref. Q, and 
by observing the oscillations of the center of the atomic 
density distribution. 



2. THE MODEL 

The quantum dynamics of an array of condensates 
in a deep enough OL, can be described by the Bose- 
Hubbard (BH) model [p. Let V H (r) = E| =1 mf2jr|/2 
be the harmonic trapping potential and let Vj,(r) = 
h 2 uj 2 sin 2 (fcri)/(4E r ) be the one-dimensional OP, where 
k is the laser wavelength, E r — h 2 k 2 /(2m) is the recoil 
energy, and lo is the angular frequency of the parabolic 
approximation of Vl at each minimum. Then, the BH 
Hamiltonian, written in terms of the boson operators a,j 
and af that, respectively, annihilate and create atoms at 
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the j site of the lattice, reads 



H = } j [ Urij (m — l)+ XjUj 



T 
2 



h.c.) , (1) 
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where the operators rii = afdi count the number of 
bosons at the i site, and the boson operators satisfy the 
standard commutation relations [cij, d^] = Si t j. The 
indices i,j € Z label the local minima rn,rij, where 
r u = iri/k, of V(r) = V L (r) + Sf =2 mOfrf/2 along the 
lattice. As the BH model describes a closed system, the 
total number of bosons N — T^jiij is a conserved quan- 
tity. Within the Gaussian approximation, the Hamilto- 
nian parameters have the following expressions in terms 
of the trapping potentials and of the optical one (see (HI). 
U := agfla^ mhQo / (2tt) is the strength of the on-site re- 
pulsion, in which we have set O = ^^^2^3- The latter 
approximation is suitable in the limit of tight confinement 
of BECs in every well. In fact, in this case the spatial 
width of each trapped condensate does not depend, in 
the first approximation, on the number of atoms in the 
well, and the condensate wave function in every potential 
minimum is well approximated by a Gaussian [15j . The 
site external potential is Xj :— e(u>) + j 2 Tr 2 h 2 ill/(4E r ), 
where e(uj) = h(ui + Q 2 + ^3)/2, and 



Hamiltonian is 



T(u) 
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1E r 



(2) 



H = Y J [U\^ + A;M 2 - ~ {Z*Z 3 + 1 + C.C. 



(4) 



where j and j + 1 run on the chain sites, and the following 
equations of motion result: 

itiZj = (2U\z J \ 2 + Xj)zj - j(zj-i + z j+1 ) , (5) 

together with the complex conjugate equations. The con- 
straint on the total number of bosons is now satisfied on 
average, the quantity N — J2j \ z j\ being conserved. 

Equations (O are the discrete version of the Gross- 
Pitaevskii equation [TH, and the corresponding super- 
fluid ground state is approximately given by the discrete 
Thomas-Fermi solution 



N 
~~M' 



(A, - A) 
2U 



(6) 



where M' — min(M,q) and q in turn is the maxi- 
mum integer such that 2UN + q(X — X q ) > 0, and 

* = 5 W „V M ' 



is the tunneling amplitude between neighboring sites. It 
is worth stressing that the Gaussian approximation is not 
essential for the use of the BH model, but it is useful in 
order to derive an analytic estimation of the Hamiltonian 
parameters. 

When the OL amplitude is weak enough so that 
T/U ^> 1, the ground-state configuration of Hamilto- 
nian (TTJ) admits a factorization into a product of site 
states that catches the superfluid nature of the sys- 
tem. The system's order-parameter dynamics near the 
ground state can be studied by a time-dependent vari- 
ational principle (TDVP) [13 ■ Following the TDVP 
method, we describe the system in terms of the trial state 
I*) = exp(iS/h)Hi\zi), which contains a product of site 
Glauber CSs. In fact, in Ref. [l6| it is shown that, in 
the limit U = 0, N, M -> 00 at fixed density N/M, the 
ground state of Hamiltonian ([1]) with M wells, is indis- 
tinguishable by a product of local coherent states. Thus, 
in the strong hopping regime, such a state should still be 
a good approximation. Here 

00 n 

\ Zi):=e -h\^Y^'^{a\T\0) (3) 

71=0 

have the defining equation a,i\zi) = and the z% 

are complex numbers. The equations of motion for 
the Zj dynamical variables are derived by a variation 
with respect to zj and z\ of the effective action S — 
J dt(iYjjijZ* — Ti) that is associated with the classical 
Hamiltonian H(Z, Z*) := (Z\H\Z), where \Z) = ^ \zj). 
Hence, the zj — (^f\aj\^) represent the classical canon- 
ical variables of the effective Hamiltonian H and satisfy 
the Poisson brackets {z*,zi} = iSje/h. The classical 



3. DYNAMICS 

After having prepared the system in the superfluid 
ground-state configuration ©, it is taken to the Mott 
regime by abruptly increasing the OP depth and, after 
an adjustable time r, it is carried back to the super- 
fluid regime by suddenly decreasing the OP depth to 
the original value. Our goal is to derive the equations 
of motion that describe the dynamics of the system af- 
ter the latter decreasing of the OP amplitude. We will 
show that these equations of motion are more compli- 
cated than the standard ones recalled in Eq. ([5]), and 
incquivalcnt to these. Furthermore, as we said above, 
the superfluid ground state is well approximated by a 
product Yit \ z i) of CSs. Hence, the condensate in each 
site i is described by a CS \zi), that is, a semiclassical 
state. On the contrary, we shall show that, after this 
stimulation, the system will be described by a product 
of integrals of site CSs. This means that the semiclassi- 
cal nature of the site states is partially destroyed during 
the intermediate quantum regime, in spite of the system 
being in a superfluid regime. 

Following the procedure above, at the time t = the 
amplitude of the OP is suddenly increased by varying 
uj from its initial value luq to ui — u\ 3> luq. Since the 
tunneling amplitude T(ui) in Eq. ([2]) is dominated by 
the exponential term exp{— 7r 2 /ia;/8£' r }, it will result in 
T(u})/T(u}q) — > for uj — ► u>\. Meanwhile, also U and Xj 
are modified when changing the OP amplitude, but their 
dependence on uj is much less dramatic. In fact, we have 
U(wi) = U y/uji/ujQ and Xj[pJx) = Xj + h(uj\ — ujq)/2. In 
order to apply the sudden approximation when the po- 
tential amplitude is varied, that is, for < t < Tb with 
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7i <C t, the jump in the potential depth must be fast com- 
pared with the tunneling time between neighboring wells, 
but slow enough so that no excitations are induced in 
each well, that is 2tt/lu, 2ir/Q, 2 , 27r/S7 3 < r b < H/T(uj ). 

For t > Tb (we will assume t\, = from now on), 
the system enters into the Mott regime and the classi- 
cal description of the system dynamics is no longer al- 
lowed; thus we resort to the quantum one. The appro- 
priate dynamics is described by the Schrodinger equation 
with Hamiltonian |T]) in which we have to set T = 0, 
U = U{u\) =: U and Aj = Xj(w\) =: Ay. The quantum 
time evolution of the initial state ([3]) is 



(7) 



where := exp[i/ft(?7 + A,-)i], = Ut/h, and the 
are those defined in ((6|). We want to stress that, al- 
though the factorization of the state vector still holds, the 
term exp[— in?u(t)] in Eq. breaks the CS structure 
of the initial state ([3]), and cij\(t)} ^ a(t)\(t)). By direct 
calculation, one can easily verify the following relations: 

((t)\aja j \(t)) = \z j \ 2 , 

Zj(t) := ((t)\aj\(t)) = Zje lht-*M a <tol*"{t)] 

e -2\ Zj \ 2 sin 2 [«(*)] ^ /gN 

((*)I4 +I %|(t)) = z ; +lZje - 2 (l^l 2 + l^ +1 | 2 ) S in 2 [«(*)] 

e t(^J-^i+i)*- ; i(|zj| 2 -|z J + i| 2 )sin[2ti (t)] _ 

The system shows the characteristic scenario of phase col- 
lapse and revivals, observed in many BEC systems [1, @. 
For < t < t, the wells' populations ((t)\af a,j\(t)) do 
not change, whereas site wave functions Zj(t) are dy- 
namically active. The modulus of Zj(t) is a periodic 
function of t, whose revival time is T m = ttH/U. The 
phase of the wave function ipj := &vg[zj(t)/\zj(t)\] = 



Xjt/h 



' sin[2w(t)] of the site j, is driven by the 



three time scales 2T m , = 2irh/Xj, and, for the sites 
where \zj\ 2 > 2ir, T Zj , the solution of the equation 
\zj\ 2 sin 2 [u(£ + T Zj )] = \ Zj \ 2 sin 2 [u(t)} +2tt. Moreover, the 
site-dependent external potentials Xj induce a dephasing 
between the wave functions of near sites: (pj — <fj+i = 
-(Trhn 1 /2) 2 (2j + l)/E r ~{\z J \ 2 ~\z J+1 \ 2 )sm 2 [u{t)}. Such 
dephasing leads the system out of the ground-state con- 
figuration. 

After a time r, the system is taken back to the super- 
fluid regime, T/U 1, by abruptly decreasing the OP 
depth (in a time of order Tb f=a 0) to its initial value. The 
t = t initial state 

K T )> = IT £ i E [ -jL=e-"^\n 3 ) (9) 



j£l M nj =0 V 

given by Eq. Q, where £j = exp{— \zj\ 2 /2}, z] 



ZjVj(r), and wi = u(t), by the identity 



lim / dx exp[— (p+e)x — inx] = exp[— n /(4p)]\/7r/p, 

with p = —i/(Aui), can be rewritten as the superposition 
of product states of CSs at each site, 



K-)> = n / 



dxj 



. e -«ir/4 e i^/(4«i)| z l e -«i) , (10) 



z^e 



j are the normal- 



where the states labeled by z'J 
ized CSs of Eq. © with Zj = z'J. " 

The time evolution for i > r of the mean-field state 
(fTU)l can be derived within the TDVP picture, in a way 
similar to that previously described, by resorting to a 
suitable superpositions of Glauber CSs. Thus we intro- 
duce the trial state = exp(iS/h)\Z) v , where \Z) V := 
Hi\zi,ui) v is written in terms of the states \zi,u\) v that 
in turn are superposition of the standard Glauber ones 
as 



\Zi,U\) v 



dxf{x,ui)\ Zi e- ix ) . (11) 



Here \z{e %x ) are the standard CSs given in ([3]), and 



f(x,m 



,-iir/i ix 2 /(4ui) 



(12) 



A remark about the trial state that we have chosen is in 
order. The TDVP method provides the best approxima- 
tion to the true state within the restricted set of states 
caught by the trial one; thus, in general, we do not know 
what superposition of the canonical CSs gives a class of 
states broad enough to obtain a good approximation of 
the true state. However, in this case, the form of the 
superposition (fTTj) is suggested by the fact that it in- 
cludes the initial condition (TIT))) . Furthermore, in the 
limit t — ► 0, that is, u\ — > 0, \zi,u\) v — + \zi) and, in 
this way, the canonical CSs ([3]) and the standard dy- 
namics ([3]) are recovered. From now on we drop the 
explicit dependence on u\ in \zi,ui) v . The scalar prod- 
uct between the \zi) v states is defined as v (zj\zi) v := 
J dxdyf*(x,ui)f(y,ui)(zje~ lx \z£e~ lv ) (here the integra- 
tions over x and y run from — oo to oo); thus, from the 
definition pip and by the normalization of the Glauber 
CSs (a\a') = exp[aW - l/2(|a| 2 + |a'| 2 )], the following 
identities can be checked by direct calculation: 

v( z j\ z j)v =1 J 



ih 

v (z j \ihd t \z j } v =—(zjZj - z*jZj) , 

v {zj\a^a]ajaj\zj) v =\zj\ A , 
v{zj\a)jaj\ Z j) v —\zj\ , 



(13) 



v { z j\ a j a j+l\ z j+l)v =ZjZj+U 



-2(|2 j | 2 + |z J + 1 | 2 )sin 2 [ii 1 ] 



J(\z j \ 2 -\z j + 1 \ 2 )sin[2u 1 ] 
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Following the TDVP procedure, we require the trial state 
to be a solution of the weaker form of the Schrodinger 
equation, 

(^\ihd t - H\y) =0, 

where dt is the time derivative and H is the BH Hamil- 
tonian (fTJ). From the latter equation, we obtain 

s = ih v (z\d t \z) v - v (z\h\z) v , 

and, from this and the relations in (I13|) . we get 



s = J dt[mJ2 - - n[z, z*)\ . (w) 

3 

The effective classical Hamiltonian H(Z,Z*) := 
V (Z\H\Z) V , can be derived by exploiting the identities 
in (fl3|) and the result is 

n{z,z*) = Y J [u\z ] \ A + \ ] \z ] \ 2 ]+ 

3 

_|^[ z * Zj+1 ^(l^l 2 -l^ +1 | 2 )-[2 Ml ] +cc] (15) 

3 

e -2(\z ] \ 2 + \z j + 1 \ 2 )sm 2 [u 1 ] _ 

The variation of the action (|14[) with respect to Zj and 
z* brings us to the classical equations of motion 

ihij = (2U\zj\ 2 + Xj)zj -yx 



i * 2 

Zj + i ZjZ j+1 



]\ Z 3? Z 3 + l 



Zj-i 



* 2 
ZaZ 



(16) 



z 2 p^+ 1 'i - -z* z 2 e^i-^ 



for j S J M , where = z*2; fe -(|z :) | 2 + |z fe | 2 )/2, and with 
the complex conjugate equations. It is worth noting that, 
in the present case, the dynamical variables Zj , z\ are 
related to the expectation values of the boson operators 
in a more complicated form than in the Glauber CS case. 
In fact we have v {zj\a,j\zj) v , = zje~~ mi exp[\zj\ 2 (e~' l2ui — 
1)]. Despite that, the \zj\ 2 still count the number of 
atoms in each site j, see the fourth of Eqs. (fie 



4. NUMERICAL SIMULATIONS 

We have numerically integrated Eqs. and (p~6|) on 
a lattice of 256 sites, with the initial conditions (fT"3f with 
the Zj given by ©, and for the values of u\ correspond- 
ing to the waiting times r = n/10, 7r/4, 7r/2. For a BEC 
in an OL with T/(2UN) as 0.01, in a harmonic trap 
with Xj/{2U) = a{j - 128) 2 and a = 0.01, the stan- 
dard equations of motion ([S]) imply a superfiuid dynam- 
ics [20(. This is shown in Fig. [1] (a), where we plot the 




FIG. 1: Center of the atomic density distribution x cm = 
Y2j \ z j\ 2 j' as a function of the rescaled time s = tfi/U (both 
dimensionless) . (a) shows the regular oscillations of the center 
of the atomic density distribution of a BEC in a harmonic trap 
obtained by numeric integration of Eqs. ((5j). (b) shows the 
same quantity obtained by numerical integration of Eqs. (|16p 
with t = 7r/10 (continuous), 7r/4 (dashed), and n/2 (dotted) 
lines. 



regular oscillations of the center of the atomic density dis- 
tribution along the chain. These oscillations have been 
triggered by displacing the condensates respect to the 
harmonic trap center, as in the experiment of Ref. [||. 
On the contrary, once the nonstandard equations 

of motion (|16|) . with the same initial conditions, entail 
insulator (dissipative) dynamics for the system. This 
is clearly shown in Fig. [T] (b), where we plot the mo- 
tion of the system's center of atomic density distribution 
for r = 7r/10 (continuous line), 7r/4 (dashed line), and 
tt/2 (dotted line). We also have performed numerical 
simulations of Eqs. (f5|) and p^|) by choosing initial con- 
ditions in the form of wave packets of Gaussian profile 
Zj = \/fccxp{— (j — x) 2 / a 2 + ip(j — x)}, where x = 128 is 
the initial center of the Gaussian, p = 3n /A is the initial 
center of mass momenta, a = 10 is the width of the Gaus- 
sian profile, and k = . exp(—(j — a;) 2 /(2cr 2 ))] -1 . The 
dynamics of these profiles have been numerically and an- 
alytically studied in Ref. [18j , where a dynamical stabil- 
ity phase diagram for these states was derived. Therein, 
and also here, the theoretical configuration where the 
harmonic trapping is off was considered, and the dynam- 
ics takes place on a finite lattice endowed with reflecting 
boundary conditions. Thus, we have performed simula- 
tions in the same conditions and we have chosen the com- 
bination of the dynamical parameters T/(2UN) — 4.17, 
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FIG. 2: Breather excitation obtained by numeric integration 
of Eqs. ([5} with a Gaussian initial profile. The breather 
travels along the chain and bounces at the lattice ends. After 
any bounce, the breather is reconstructed. 
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i 

FIG. 3: Excitation obtained by numeric integration of Eqs. 
(|16l) with the same initial condition as in Fig. [5] The value 
of the waiting time is r = n/2. This excitation is completely 
destroyed after a few bounces at the lattice ends. Thus, in the 
far-from-equilibrium situation, this excitation loses stability 
and it seems to behave like the states of the diffusion regime 
that has been identified in [l8| ]. 

which corresponds to the region of the phase diagram of 
Ref. [IH where (stable) breather excitations were identi- 
fied. In Fig.[5]we report the two-dimensional contour plot 
obtained by numeric integration of ([5]) with this Gaus- 
sian initial condition. In Fig. [2] is clear that the breather 
structure is maintained when the traveling excitation is 
reflected at the lattice boundaries. Figure [3] shows the 
same quantity obtained by integration of Eqs. (|16p with 
the same initial condition as in Fig. [21 and with r = 0.5. 
Figure [3] clearly shows that the initial excitation, inte- 



grated with nonstandard dynamics, pretty soon loses sta- 
bility, emitting atoms incoherently at any bounce with 
the lattice boundary. 

5. FINAL REMARKS 

In the present paper we have studied how the dynam- 
ics of a superfluid is affected by briefly bringing the sys- 
tem into the insulating regime. We have shown that the 
system is taken to an excited state, described by a super- 
position of product states of Glauber coherent states at 
each site, which we have derived. The classical equations 
of motion ruling its dynamics have been derived. Fur- 
thermore, we have shown that these classical equations 
of motion are inequivalent to the standard discrete non- 
linear Schrodinger equations that describe the dynamics 
of an array of BECs [y] in the superfluid regime. By nu- 
merically integrating such nonstandard equations with 
several initial conditions, we have shown that the system 
loses coherence, becoming insulating. 

The simulations we have performed show that the in- 
terplay between classical and quantum dynamics leads 
to loss of the coherence properties of the system. In 
fact, the brief period in the insulating regimes changes 
the superfluid wave function, which becomes a superpos- 
tion of product states of the site's coherent states, that 
is, a product of mean-field states. Each mean-field state 
has a complicated distribution of phases at each site that 
results from the intermediate quantum dynamics. This 
distribution of phases leads to a unique tunneling dy- 
namics described by a complicated hopping term. By a 
glance at Eqs. (fTl)]) one can guess that, as a consequence 
of this unusual term, the "effective" tunneling rate be- 
tween close sites in the case of Eqs. (fT6)l becomes site 
(population) and time dependent. For this reason the 
system loses coherence. 

It is also worth emphasizing that the nonadiabatic 
procedure we have described in the present paper can 
straightforwardly be realized in a real experiment similar 
to those of Refs. [H, [Hj]. Therefore, by displacing the 
contensates with respect to the harmonic trap, as done 
in the experiment of Ref. || , and by observing the oscil- 
lations of the center of the atomic density distribution, 
the effects of the nonstandard dynamics can be directly 
observed. 
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